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Abstract. This paper proposes a simple and efficient method for the 
reconstruction and extraction of geometric parameters from 3D tubu¬ 
lar objects. Our method constructs an image that accumulates surface 
normal information, then peaks within this image are located by track¬ 
ing. Finally, the positions of these are optimized to lie precisely on the 
tubular shape centerline. This method is very versatile, and is able to 
process various input data types like full or partial mesh acquired from 
3D laser scans, 3D height map or discrete volumetric images. The pro¬ 
posed algorithm is simple to implement, contains few parameters and 
can be computed in linear time with respect to the number of surface 
faces. Since the extracted tube centerline is accurate, we are able to de¬ 
compose the tube into rectilinear parts and torus-like parts. This is done 
with a new linear time 3D torus detection algorithm, which follows the 
same principle of a previous work on 2D arc circle recognition. Detailed 
experiments show the versatility, accuracy and robustness of our new 
method. 


1 Introduction 

Tubular shapes appear in various image application domains. They are com¬ 
mon in the medical imaging field. For instance, blood vessel identification and 
measurements are an important object of study [nnzEa. The wall thickness in 
bronchial tree plays also an important role in several lung diseases m- Tubular 
shapes also occur in CT volumetric images of wood m- their segmentation into 
knots is exploited by agronomic researchers or in industrial sawmills. Outside 
volumetric images, tubular objects are also present in industrial context with the 
production of metallic pipes from bending machines. Quality assessment of such 
metallic pieces is generally achieved with a direct inspection by a laser scanner. 
Such process is also performed for calibration purpose and reverse engineering 
tasks. 

* This work was partially supported by the ANR grants DigitalSnow ANR-11-BS02- 
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Geometric properties of tubular structures are extracted in different ways 
depending on the application domain and on the nature of input data. From 
unorganized set of points, Lee proposed a curve reconstruction exploiting an 
Euclidean minimum spanning tree with a thinning algorithm and applies it to 
pipe surface reconstruction m- Later, Kim and Lee proposed another method 
based on shrinking and moving least-squares m to improve the reconstruction 
of pipes with non constant radius. However, as shown by Bauer and Polthier 
[6], such a reconstruction method produced noisy curves in particular for data 
extracted from partial scans like the ones of Fig.[^ (b). Another approach esti¬ 
mates the principal curvatures of the set of points in order to detect cylindrical 
and toric parts m- Although promising, this approach suffers from the quality 
of the local curvature estimator. To overcome this limitation, Bauer and Polth¬ 
ier [6] proposed to recover a parametric model based on a tubular spine. Their 
method is able to process partial laser scans limited to one particular direction. 
The main steps of their method consists in first projecting the mesh points onto 
the spinal region of the mesh before reconstructing a spine curve and analyz¬ 
ing it. The method requires as parameter one radius size, and it cannot process 
volumetric data (voxel sets) or heightmap data. 



Fig. 1. Different kinds of 3D tubular data: (a) input data obtained from full laser scan, 
(b) partial scans from one direction, (c) digital set of voxels and (d) height map. 


More generally, classic medial axis extraction looks to be a potential solu¬ 
tion for tubular shape analysis [S]. However such extraction may be sensitive to 
noise or to the presence of small defaults in the volumetric discrete object (like 
small hole). Fig.shows some results obtained with different methods available 
from the implementation given in the authors survey. We can clearly see that 
small holes in the digital object significantly degrade the result. More recently, 
many advances came from the field of mesh processing with approaches based 
on mesh contraction However, they are generally not adapted to surface 

with boundaries like the partial scan data of Fig.[2 In the same way, they are 
not simple to adapt to volumetric data like digital object made of voxels or 
height map. To process volumetric discrete objects. Gradient Vector Flow m 
was exploited by Hassouna and Farag in order to propose a robust skeleton curve 
extraction m- This method was also adapted to process gray values volumetric 
images used in virtual endoscopy [5]. In the field of discrete geometry we can 
mention a method which propose to specifically exploit 3D discrete tools to ex¬ 
tract some medial axis on grey-level images [ 3 ] . To sum up approaches on medial 
axis, they are designed to process shapes defined as volumes, but they fail when 
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processing open surfaces or partial samplings of the shape boundary. To process 
such data, we can refer to the work of Tagliasacchi et al. [3T], who propose an 
algorithm based on surface normals. However the resulting quality depends on 
manual parameter tuning. 

In this work, we propose a unified approach to the reconstruction and the ge¬ 
ometric analysis of tubular objects obtained from various input data types: laser 
scans sampling the shape boundary with partial or complete data (Fig. [I](a,b)), 
voxel sets sampling the shape (Fig.[^(c)) or more specifically from height map 
data (Fig.[^(c,d)). Potential applications of the latter datatype are numerous be¬ 
cause of the increasing development of Kinect@-\ike devices. Our main contribu¬ 
tions are first to propose a simple and automatic centerline extraction algorithm, 
which mainly relies on a surface normal accumulation image. Like other Hough 
transform based applications um, this algorithm can handle various types of 
input data. We also propose to extract geometric information along the tubular 
object, by segmenting it into rectilinear and toric parts. This is achieved with a 
3D extension of a previous work on circular arc detection along 2D curves. In 
the following sections, we first introduce the new method of centerline detection, 
then we show how to reconstruct the tubular shape and decompose it into mean¬ 
ingful parts. We conclude with representative experiments showing the qualities 
of our method. 



(a) thinning 
1 min 4s 


(b) geometric (c) potential field 
0.5 s 3min 


(d) proposed 
6s+3s 


Fig. 2. Skeleton extraction from three different methods presented in [8] with the 
implementation given by the authors. Our method is presented on the right. 


2 Fast and Simple Centerline Extraction on 3D Tubular 
Object 

In this section, we present centerline extraction algorithm based on surface nor¬ 
mal accumulation. It consists in three main steps. First, we compute a 3D accu¬ 
mulation image, which counts for each voxel how many faces of input data have 
their normal vector pointing through the voxel (i). Depending on input data, 
normal vectors can be defined directly from the mesh faces or estimated by a 
more robust and accurate estimator (in particular if we process a digital object). 
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Then a tracking algorithm (ii) extracts an approximate centerline by following 
local maxima in the accumulation image. Finally, to remove digitization effects 
due to the 3D discrete accumulation image, we optimize the position of center 
points (iii) by a gradient descent method. 

2.1 Accumulation Images Prom Normal Vectors 

The first algorithm requires as input a set of faces (a mesh or a digital surface), 
their associated normal vectors, and a 3D digital space (the 3D grid that will 
store the accumulated values). If the input object is a mesh, the gridstep of 
the digitization grid must be specified by the user. Depending on the awaited 
accuracy, a default gridstep can be chosen as the median size of mesh faces (with 
the face size defined as its longest edge). If the input object is a digital object 
or a heightmap, the digitization grid just matches their resolution. Besides, this 
algorithm requires as parameter some approximation of the tube radius R. 

The whole algorithm is detailed in Algorithm[^ It outputs for each voxel the 
number of normal vectors going through it as well as a vector estimating the 
tube local main directions (i.e. the tangent to the centerline or equivalently the 
direction of minimal curvature along the tube boundary). Fig.illustrates the 
main steps of the algorithm with the 3D directional scans, starting from the 
face origin fk in the direction of its normal vector along a distance denoted 
accRadius (set to i? + e where e is used to take into account possible small 
variations of the radius along the tube, see image (a) of Fig.Q. During the 
scan, the accumulation scores are stored for each visited voxel (image (b) of 
the same figure). The principal direction of a voxel is also updated for each 
scan (image (c)). More precisely, if we denote by hy and Uk the two last normal 
vectors intersecting a voxel V for the scans j and k, the principal direction dk 
for the current scan is given by: dk = dj + (h^ order to ignore non 

significant directions induced by near colinear vectors, we add a small constant 
(set by default to 0.1) to filter the norm of the resulting vector A 



(a) directional scan (b) accimage 


(c) dirlmage (d) lastVectors 


Fig. 3. Illustration of Algorithm[^ which builds an accumulation image whose peaks 
match the centerline of the tubular shape. 


Fig.|4] illustrates the computations made in Algorithmic for a mesh input 
data, and it shows the resulting accumulation image (image (c)) and vectors 
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it orthogonal to the tube main direction a . Since centerline extraction relies 
on these accumulation images, we evaluate the robustness of these images with 
various input surface types. The first row of Fig.|^ presents the resulting 3D 
accumulation images obtained on partial, noisy, digital or on small resolution 
mesh. In all these configurations, relative maximal values are indeed well located 
near the center of the tubular shape. A fixed threshold was applied in order to 
highlight the voxels with accumulation values close to maximal ones. Such voxels 
are drawn in black and for a particular selected voxel, we have highlighted their 
scanning origin faces with blue lines. All these results confirm the robustness of 
the proposed algorithm. We have therefore a solid basis for the tracking algorithm 
presented in the following part. 



Fig. 4. Illustration of accumulation images generated from surface normal vectors. 
Image (a) (resp. (b)) illustrates some (resp. all) scanning directions defined from mesh 
triangles. Image (c) illustrates the values obtained in the 3D accumulation image and 
(d) shows some voxels having accumulation score upper than a given threshold. In the 
last image we also display the set of faces contributing to the accumulation score. 


2.2 Centerline Tracking from Image Accumulation 

Even if the maximal values obtained in the pre¬ 
vious part are well centered on the tubular object, 
a simple thresholding is not robust enough to ex¬ 
tract directly the centerline. Furthermore it implies 
the manual adjustment of the threshold parameter. To illustrate this point, the 
image on the side shows different results obtained by choosing various threshold 
parameters a. A too strict threshold implies disconnected points, while a less 
restrictive one produces a thick line with parasite voxels. 

To better approach the centerline we propose to define a simple tracking al¬ 
gorithm exploiting the output of Algorithmic i.e. the accumulation image and 
the direction vectors image. As described in Algorithmic the main idea is to 
start from a point Cq detected as a maximal accumulation value of the 3D 
accumulation image. Then, from a current point Ci of the centerline, the algo¬ 
rithm determines next point as the point having maximal accumulation 

value in the 2D patch image defined in the plane normal to the direction 

dirlmage((7i) at distance trackStep (see Fig. |6](a,b)). 
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Fig. 5. Experiment of the robustness of the surface normal accumulation algorithm 
applied on different types of input surface: on sector filtered mesh (a,f), on partial scan 
mesh (b,g), on noisy mesh (c,h), on digital surface (d,i), and on small resolution mesh 

(ej). 


Algorithm 1: accumulationFromNormalVectors : From position and nor¬ 
mals of faces of an input mesh, this algorithm computes an accumulation 
image (accimage) by a directional scan starting from a face center in the 
direction of the face inward normal. It also outputs the image of vectors 
representing the local main axis direction of the tubular shape (dirImage). 

Input : mesh // Triangular mesh of a tube 

accRadius // Accumulation length from center of faces 
minNorm = 0.1 // Minimum norm value 

Ouput : accimage / / Accumulation of normal vector number passing through a coordinate 
dirlmage // Cross product of all normals passing through a coordinate 
maxAcc // Maximum number of normals passing through an (x,y,z) coordinate 
maxPt // maxAcc coordinates 

Variable: lastVectors // The last considered normal for each (x,y,z) coordinate 

mainAxis // Vector contributing to the cross product of a directional vector 
lastVectors = Image3D(mesh.dimensions()) 
maxAcc = 0 

foreach face in mesh do 

currentPt = face.center 

normalVector = face.normalVector().normalizedO 
while distance (currentPt, face, center) < accRadius do 
if accimage [currentPt] != 0 then 

mainAxis = lastVectors[currentPt] X normalVector 
if norm (mainAxis) > minNorm then 

|_ dirlmage[currentPt] += mainAxis*sign(mainAxis • dirlmage[currentPt]) 

lastVectors[currentPt] = normalVector 
accimage[currentPt]++ 
if accimage [currentPt] > maxAcc then 
I maxAcc = accimage[currentPt] 

L maxPt = currentPt 
L currentPt += normalVector 
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Fig. 6. Tracking algorithm step. The patch generated from the maximum Ci 

of I patch fhe dk direction at a trackStep distance before to localize the maximum 

Ci+i of this new patch. 


Algorithm 2: trackPatchCenter: tracking algorithm in one direction, 
given a starting point and an orientation. 

Input : accimage // Accumulation of normal vector number passing through a coordinate 
dirlmage // Cross product of all normals passing through a coordinate 
accRadius // Accumulation length from center of faces 
startPt, // Start point for tracking (must belong to the centerline) 
trackInFront // True if tracking direction is in the startPt vector direction 
trackStep // Distance between two consecutive centerline points 
Output : centerline // Point set constituting the tube centerline 
Variable: continueTracking // True if tracking can continue 
patchSize // Dimension of the square patch 

currentPt, previousPt // Considered point during an iteration 
lastVect // Directional vector associated to previousPt 
centerPatch // Patch center finding from currentPt 
patchImageSize = 2 * aRadius ; 
centerline = emptySetO 
continueTracking = true 
patchSize = 2 ♦ accRadius 
currentPt = startPt 

lastVect = trackInFront ? dirlmage( startPt ) : - dirlmage( startPt ) 
previousPt = startPt - lastVect ♦ trackStep 
while continueTracking do 

centerline.append( currentPt ) 
dirVect = dirlmage[currentPt].normalizedO 
if lastVect .dot (dirVect) < 0 then 
L dirVect = -dirVect 

continueTracking = isInsideTube( accimage, currentPt, previousPt, trackStep, 7r/3 ) 
previousPt = currentPt 

// Defined the next image patch center point 
centerPatch = currentPt + ( dirVect * trackStep ) 
if not accimage. domainO . contains ( centerPatch ) then 
L break 

// Extract a 2D image of size 2 * accRadius from the 3D image accimage, centered on 
// centerPatch and directed along dirVect 

patchimage = extractPatch( accimage, centerPatch, dirVect, 2 ♦ accRadius ) 
maxCoords = getMaxCoords( patchimage ) 
lastVect = dirVect 
previousPt = currentPt 

_ currentPt = patchSpaceToAccImageSpace( maxCoords) 
return centerline 
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2.3 Skeleton position optimization 

Since the resulting tracking skeleton is embedded in a digital space, it suffers 
from digitization artefacts and is not perfectly centered within the input mesh. 
Moreover, depending on normal mesh quality, the tracking algorithm can poten¬ 
tially be influenced by perturbated normal directions, and may deviate from the 
expected centerline. Such perturbations can dramatically degrade the quality of 
upcoming geometric analysis, and hence impose some unwanted post processing 
tasks. To avoid such a difficulty, we propose to apply an optimization algorithm 
in order to obtain a perfectly centered spine line. 

The idea is to model the quality of the current fitting by an error Es{C), 
defined as the sum of the squared difference between the known tube radius R 
and the distance between the tube center C and its associated input mesh points 
Mi. We wish to find the best position for center C that minimizes this error. 
Otherwise said, we look for the circle of radius R that best fits the data points 
Mi in the least-square sense. Hence, the error is 


AT-l 

Es{C)=J2{\\^i\\-Rf. ( 1 ) 

i=0 

This minimization problem is easily solved by a gradient descent algorithm that 
follows the direction of steepest descent of the error. By simple derivation, its 
gradient is 

N-l 

The gradient descent can also be interpreted as elastic forces acting on the center 
C and pulling or pushing it in the direction of data according to the current 
distance. Then the minimization process applies at each step of the process the 
sum of theses forces on (7, giving with the notations of Fig. 07 

By this way, at each step, the total error Es decrease and we iterate the 
process until convergence, i.e. the difference of errors between two iterations is 
below a fixed e^. 

Contrary to a simple average of neighborhood points, the optimization per¬ 
forms well even on partial mesh data, with missing parts or holes. Moreover, it is 
possible to ponder each force with its face area in order to better balance forces 
in presence of irregular sampling with variable density. 

3 Reconstruction Results and Geometric Analysis 

Reconstruction results. Several centerline extractions and tubular shape re¬ 
constructions are shown on Fig.|^ Our input dataset contains several types of 
metallic tubes numerized with different acquisition tools. In each case, the cen¬ 
terline is always well delineated without the need to tune a special parameter. 
When the input data is a complete or partial mesh, the normal is simply esti¬ 
mated as the cross product of face edges. When the input data is a digital object 
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Fig. 7. (left) Illustration of the process to optimize the centerline position (point C) 
with elastic forces (blue arrows). Each elastic force is attached to one point of the 
input mesh sector (a point Mi represented in black) and oriented in the direction of 
the center of the virtual circle of center S. (right) evolution of the convergence speed 
in the optimization process. 


or a height map, we use the digital Voronoi Covariance Measure [9] to estimate 
the normal vector. Parameters of this estimator are easily set since we know the 
radius of the tubular object; they also have little influence on the result, since the 
accumulation image makes the process very robust. The running time was less 
than 30s for each experiment. As show on image (j) and (f), few places present 
reconstruction errors. Small errors may be found near bent areas. These errors 
are less related to the reconstruction method than to the physical shape under 
study, since the bending machine has deformed the tube at these places. 


Geometric analysis with 3D 
arc detection. We wish to seg¬ 
ment the tubular shape into rec¬ 
tilinear and toric parts. This pro¬ 
blem is equivalent to segment- (a) input polygonal curve representa- 

ing its centerline into 3D straight 

segments and 3D circular arcs. We thus extend to 3D a method presented by 
Nguyen et al. [18], which was designed to cut a 2D discrete curve into straight 
and circular pieces. It relies on properties of circular arcs in the tangent space 
representation that are inspired from Arkin [2| and Latecki M- The tangent 
space representation of a sequence of points C = defined as follows : 

Let li be the length of segment CiCi^i and ai = Z(Ci_iCi, CiQ+i). Let us 
consider the transformation that associates C with a polygon of constituted 
by segments To 2 T ('.-1 i^i, 0 < i < n (upper floating figure) with: 

rL = (0,0), Ta = (’ri.i:;.!;"om i to „Jt., = 

{Tii.x^Tii.y + ai) for i from 1 to n — 1. Moreover, let M = (M^)^^q^ be the 
sequence of midpoints of segments for i from 0 to n — 1. The main 

idea of the arc detection method is that if (7 is a polygon that approximates a 
circle or an arc of circle then {Mi)^~Q is a sequence of (approximately) co-linear 
points m- 
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(a) 

(b) 

(c) 



Full mesh 

R = 6, t=6.23s, 151 444 faces 



(g) 



(k) 

(l) 

(m) 

(n) 



Reduced scaned area 
R = 6, t=4.36s, 52 914 faces 


(O) 



Digital object 

R = 4.9, t=6.91s, 60 768 faces 



Height map 

R = 8,t=22.33s, 645 450 faces 


(q) straight (in light blue) and toric 
(in dark blue) segments. 

R = 6, t=12.17s, 187 638 faces 


Fig. 8. Result of reconstruction from various input data types. Images (a,g,l,p) show 
the centerline by transparency through the input surface, (k) and (o) are some input 
surfaces. (b,e,i,m) images are the reconstructed tubes built from the centerlines. Images 
(j) and (f) show the local squared distance error between the reconstructed tube and the 
input shape. Image (q) illustrates the decomposition of a tube into rectilinear and toric 
parts. For all experiments, the tracking parameter and epsilon were set respectively to 
R and 0.001. Running times correspond to executions on a MacBook computer with a 
2,5 GHz Intel Core il processor. 
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In our work, the sequence (Ci)i of skeleton points obtained in Section 2.3 


is considered and the representation of this sequence of points in the tangent 
space is computed. If the angles for k from p to of consecutive points of 
{Ck)\^p are close to 0, these points belong to a straight line. Otherwise, the co- 
linearity of the corresponding midpoints is tested in the tangent space 


iO UCDUCU. ill UiiC UCtiigCliU 

:.[^ (q) shows an example of tubular 


by using an algorithm presented in HH]. Fig.|_ 
shape decomposition with this 3D variant of circular arc detection. Toric and 
rectilinear parts are correctly identified. 


4 Conclusion and Discussion 

A new efficient and simple method was presented to solve the problem of de¬ 
lineating the centerline of 3D tubular shapes, for various types of input data 
approximating its boundary: mesh, set of voxels or height map. The method is 
robust to missing parts in input data as well as perturbations: in these situations, 
it still returns accurately the centerline position. To achieve this, we have de¬ 
composed the process in three steps : 1) computing an accumulation map from 
faces and their normal vectors, 2) tracking of centerline through cross-section 
maximas of the accumulation map and 3) optimization of the centerline position 
by a better fitting of the model to the nearest faces along the centerline. The cen¬ 
terline was accurate enough to allow further geometric analysis. We have shown 
how to decompose the tubular shape into rectilinear and toric parts by a sim¬ 
ple adaptation of a 2D circular arc detection algorithm. The hypothesis about 
constant radius parameter R only influences the skeleton position optimization. 
This limitation could be resolved either by direct radius estimation from the ac¬ 
cumulation image or by radius optimization during position optimization. This 
is left for future works. The whole process was implemented with the DGtal [T] 
framework and will soon be available in its companion DGtalTools. 
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